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This paper presents a quasi-local method of studying the physics of dynamical black holes in nu- 
merical simulations. This is done within the dynamical horizon framework, which extends the earlier 
work on isolated horizons to time-dependent situations. In particular: (i) We locate various kinds of 
marginal surfaces and study their time evolution. An important ingredient is the calculation of the 
signature of the horizon, which can be either spacelike, timelike, or null, (ii) We generalize the calcula- 
tion of the black hole mass and angular momentum, which were previously defined for axisymmetric 
isolated horizons to dynamical situations, (iii) We calculate the source multipole moments of the black 
hole which can be used to verify that the black hole settles down to a Kerr solution, (iv) We also study 
the fluxes of energy crossing the horizon, which describes how a black hole grows as it accretes matter 
and /or radiation. 

We describe our numerical implementation of these concepts and apply them to three specific test 
cases, namely, the axisymmetric head-on collision of two black holes, the axisymmetric collapse of a 
neutron star, and a non-axisymmetric black hole collision with non-zero initial orbital angular mo- 
mentum. 
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I. INTRODUCTION 

In spite of fundamental advances in our understand- 
ing of black holes, relatively little is known about them 
in the fully non-perturbative, dynamical regime of gen- 
eral relativity. Most of our intuition regarding black 
holes comes from studying the stationary axisymmet- 
ric Kerr-Newman solutions, and perturbations thereof. 
This, along with post-Newtonian calculations which 
treat the black hole as a point particle, are usually ade- 
quate for understanding many astrophysical processes 
involving black holes. However, understanding the 
gravitational waveforms arising due to, say the merger 
phase of the coalescence of two black holes or the grav- 
itational collapse of a star, will require us to go beyond 
perturbation theory and to confront the non-linearities 
and dynamics of the full Einstein equations. This regime 
may contain qualitatively new, non-perturbative fea- 
tures. In this paper, we discuss an important ingredi- 
ent for understanding this regime, namely, the dynam- 
ics of the black hole horizon. Numerical simulations 
of black holes have greatly improved in the last few 
years. Simulations of the entire merger process, start- 
ing from the last few orbits of the inspiral right up to 
the ringdown have become possible in the past year 
dSy-H.la.lalElal- It is then important to look for better 
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ways to extract more physical information from simu- 
lations and to compare results from two different sim- 
ulations performed using different coordinate systems, 
gauge conditions etc. This can be a non-trivial task in 
itself, and understanding black hole horizons is a neces- 
sary ingredient. 

Due to their global nature, black hole event hori- 
zons can only be located once a simulation is complete 
and we have obtained the full spacetime. In numeri- 
cal simulations, it is instead common to use marginally 
trapped surfaces to locate black holes on a Cauchy sur- 
face in real time. We use the formalism of dynami- 
cal horizons |9j, flOIT to study black holes. Using iso- 
lated/dynamical horizons, it is shown that marginally 
trapped surfaces, while not a substitute for event hori- 
zons, do have many useful properties and can be used 
fruitfully to study black hole physics. Dynamical hori- 
zons are a significant extension of the isolated horizon 
framework jlj |l3 Bill El, which models isolated 
stationary black holes in an otherwise dynamical space- 
time. Both these frameworks are, in turn, very closely 
related to and motivated by the earlier work on trapping 
horizons by Hayward El. EE Ell See EHHHj for 
reviews. Information obtained from these quasi-local 
horizons complements the information obtained from 
the event horizon. Once a simulation is complete and 
ready for post-processing, event horizons are useful for 
studying global properties and the causal structure of 
the spacetime, and also phenomena such as the topol- 
ogy change of the horizon during a black hole coales- 
cence. Reliable and computationally efficient codes are 
now available for locating event horizons (see e.g. 1 22]). 
Such information cannot be obtained at the quasi-local 
level, which is instead better for tracking the physical 



parameters and geometry of a black hole in real time. 

The dynamics of apparent and event horizons have 
been numerically studied in the pas t in detail in axisym- 
metry (see e.g. 1 23, HIHIMII^S]). We want to extend 
this work to non-axisymmetric and non-vacuum space- 
times, and we want to emphasise non-gauge-dependent 
analysis methods. In particular, we consider the follow- 
ing applications: (i) We study the behavior of various 
marginally trapped surfaces under time evolution. This 
leads to greater insights about the trapped region of a 
spacetime. An important ingredient here is the signa- 
ture of the world tube of marginally trapped surfaces. 
This world tube is known to be null for isolated hori- 
zons, and more generally, it can be either spacelike or 
timelike; we show that both types occur frequently in 
numerical simulations, (ii) We give meaningful defini- 
tions for the angular momentum, mass, and higher mul- 
tipole moments for the dynamical black hole. The mul- 
tipole moments capture gauge invariant geometrical in- 
formation regarding the horizon geometry, and should 
be useful for understanding fundamental issues such as 
the final state of black hole collapse. For example, we 
would expect that after a black hole has formed and set- 
tled down, its multipole moments should be identical to 
the source multipoles of a Kerr black hole. We show that 
it is, in principle, possible to verify this conjecture and to 
calculate the rate at which a black hole approaches equi- 
librium, (iii) We also describe and implement methods 
for calculating the energy flux falling into the horizon. 
This gives us detailed information on how black holes 
grow as they swallow matter and radiation. 

This paper is organized as follows. Section [n] sets 
up notation, and summarizes the basic definitions and 
properties of trapped surfaces and dynamical horizons. 
Section|ITT]describes the various physical quantities that 
we calculate using dynamical horizons, and also their 
numerical implementation. Section IIVI presents three 
concrete, well known numerical examples where these 
concepts are applied and finally, section |V| discusses 
some open issues and directions for further work. Un- 
less mentioned otherwise, we use geometrical units 
with G = c = 1, the spacetime signature is (—,+,+,+), 
all manifolds and fields are assumed to be smooth, and 
the Penrose abstract index notation is used throughout. 
The derivative operator compatible with the spacetime 
metric g a y is V fl and, following Wald |29], the Riemann 
tensor is defined via (V fl Vj, — Vf,V fl )o; c = R^/cv^. 



II. BASIC NOTIONS AND DEFINITIONS 

A. Trapped surfaces and apparent horizons 

Let S be a closed, orientable spacelike 2-surface in 
a 4-dimensional spacetime {M,g a },). The expansion of 
any such surface can be defined invariantly without any 
reference to a time slicing of the spacetime. Since S is 
smooth, spacelike, and 2-dimensional, the set of vec- 



tors orthogonal to it at any point form a 2-dimensional 
Minkowskian vector space. Thus, we can define two lin- 
early independent, future-directed, null vectors t a and 
n" orthogonal to S such that 
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(2.1) 



Note that this convention is different from that used in 
ilflH . We shall assume that we know a 'priori what the 
outgoing and ingoing directions on M. are. By conven- 
tion, i a will denote an outgoing null normal and n a an 
ingoing one. The null normals are specified only up to a 
boost transformation 



t -> ft , 



/ 



V 



(2.2) 



where / is a, positive definite, smooth function on S. All 
physical quantities must be invariant under this gauge 
transformation. 

The Riemannian 2-metric q a \, on S induced by the 
spacetime metric g a \, is 



Rab = gab + £an b + n a £ b . 



(2.3) 



The tensor q a can be viewed as a projection operator on 
to S. The null expansions are 
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(2.4) 



These expansions tell us how the area element of S 
changes as it is deformed along £ a and n a respectively. 

The shear of £", cr^\ ab , is the symmetric trace-free part 
of the projection of V a ^b : 



■c xdx 



C{l>)ab - faHb^(c^d) 

Similarly the shear of n a is 
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(2.5) 



(2.6) 



Note that these definitions only involve derivatives tan- 
gential to S. Thus £" and n" can, if necessary, be ex- 
tended arbitrarily away from S while computing these 
quantities. 

The closed 2-surface S is said to be a trapped surface 
if both expansions ©(« and &( n \ are strictly negative. 
This is very different from a sphere in normal flat space 
which has positive outgoing expansion and negative in- 
going expansion. This definition was first introduced 
by Penrose [30], who recognized its importance in the 
formation of singularities. On a marginal surface, one 
of the two null expansions vanish. Of particular in- 
terest are the marginally outer trapped surfaces (MOTSs), 
for which the outgoing null rays along £ a have zero ex- 
pansion. In addition, we shall mostly deal with future 
marginally outer trapped surfaces (FMOTSs), i.e., MOTSs 
with ©(„) < 0. 

There are three main reasons why closed trapped sur- 
faces are important for studying black holes. First, the 



existence of a trapped surface implies the existence of 
a singularity in the future 130113 ill . Secondly, they are 
guaranteed to always lie within the event horizon. Fi- 
nally in stationary spacetimes, the null generators of the 
event horizon have zero expansion. Thus for stationary 
spacetimes, the cross-section of the event horizon is a 
MOTS. 

While trapped and marginally outer trapped surfaces 
are defined in the full four dimensional spacetime, in 
numerical relativity, one usually considers trapped sur- 
faces in conjunction with a foliation of (partial) Cauchy 
surfaces containing S; it is numerically much easier to 
look for closed surfaces on the Cauchy surface rather 
than in the full spacetime manifold. For concreteness, 
we shall work in the ADM formalism where the rele- 
vant portion of spacetime is foliated by spacelike sur- 
faces, and E shall denote one of the leaves of this folia- 
tion. However, it will be obvious that the formalism is 
applicable no matter how Einstein's equations are im- 
plemented. 

The trapped region 7^ on E is defined to be the set of 
points in E through which there passes a trapped sur- 
face contained entirely in E. Note that there could be 
points in E not contained in 7^, but through which there 
passes a trapped surface not contained in E. Thus, 7^ is 
a subset of the intersection of E with the 4-dimensional 
trapped region in the full spacetime. A connected com- 
ponent of the boundary of 7^ is called an apparent hori- 
zon (AH). Under suitable regularity conditions, the AH 
can be shown to be a MOTS |32, 33]. Thus, an appar- 
ent horizon is the outermost MOTS on E. Due to this 
"outermost" property, an AH is not a quasi-local object 
on E. The behavior of AHs under time evolution can be 
quite irregular. For example, they can "jump" discontin- 
uously. On the other hand, as we shall soon see, MOTSs 
are more regular. 



B. Dynamical horizons 

1. Definition and examples 

We can use marginal surfaces to extract physically in- 
teresting information about the black hole. The key idea 
is to look not at a single MOTS by itself, but rather a 
world tube H of MOTSs constructed by stacking up the 
MOTSs obtained by time evolution. Such a world tube is 
called a Marginally Trapped Tube (MTT). An MTT is thus 
a smooth 3-surface foliated by MOTSs. 

The existence of MTTs: Numerically, it has been ob- 
served that marginal surfaces (though not apparent 
horizons — see below) usually behave smoothly under 
time evolution and produce a smooth MTT. This obser- 
vation is placed on a more rigorous footing by the re- 
cent result of Andersson et al. 1 34], which proves the lo- 
cal existence of MTTs for a large class of MOTSs. Their 
results require the MOTS to be strictly-stably-outermost. 
An MOTS S on E is said to be strictly-stably-outermost 



if there exists an infinitesimal first order outward defor- 
mation which makes S strictly untrapped. Working with 
a radial coordinate r on E such that S is a level set of 
r, and r increases in the outward direction, a sufficient 
(but not necessary) condition for S to be strictly-stably- 
outermost is d r &i(\(r) > everywhere 1 on S. Here it 
is understood that we obtain 0/« as a function of r by 
calculating &n\ for the constant-r surfaces in the vicin- 
ity of S. In principle, for an unfortunate choice of r, it 
might happen that dr®u) < even though there is a 
different choice for which this condition is satisfied. In 
any case, this is sufficient for verifying that S is strictly- 
stably-outermost. This condition, unlike the outermost 
condition for an AH, is a quasi-local condition. We have 
found in our simulations that most physically interest- 
ing MOTSs, such as ones which asymptote to the event 
horizon, and also AHs, satisfy this condition quite gen- 
erally. However, as we shall see, there exist also MOTSs 
which are not strictly-stably-outermost. In practice, in- 
stead of checking 3 r Om > directly, we look for a sur- 
face with a small positive (or negative) non-vanishing 
expansion, and check that it lies completely outside (or 
inside) the MOTS. 

It is shown in ] 34] that if a MOTS S is strictly-stably- 
outermost, then at least locally in time, S is a cross- 
section of a smooth MTT. More explicitly, this result 
shows that given a foliation of the spacetime by Cauchy 
surfaces Ej, if there is a MOTS So on Eg which is strictly- 
stably-outermost, then MOTSs St exist on Ef for — e < 
t < e (for sufficiently small e) such that the union (J St 
is a smooth MTT. The MTT will exist for at least as long 
as the MOTS remains strictly-stably-outermost. This is 
a conceptually important result for numerical relativity 
because it shows that a large class of MOTSs behave reg- 
ularly under time evolution. How is this to be reconciled 
with the known fact that AHs can "jump" during a time 
evolution? The reason is simply because of the outer- 
most property. It is possible that a new MOTS can ap- 
pear on the outside of a given MOTS. The "old" MOTS is 
then no longer the globally outermost one even though 
it is locally outermost, and it continues to evolve in a 
perfectly regular manner, but it is no longer an AH. 

There are, as yet, no similar existence proofs for 
MOTSs which are not strictly-stably-outermost. How- 
ever, as we shall see later, we find in all the examples 
we have looked at, that MOTSs evolve smoothly even 
in this case, forming a regular world tube. 

Isolated and dynamical horizons: An MTT is null in equi- 
librium situations when no matter or radiation is falling 
into it; the rest of the spacetime is still allowed to be 



1 More precisely, d r ®m{r) > with 9 r ©m(r) > somewhere on S. 

2 It is harder to show that a MOTS is not strictly-stably-outermost. 
This can be done by calculating the signature of the horizon (see 
below) or by calculating the principle eigenvalue of the stability op- 
erator defined in [34]. 



highly dynamical. This situation is formalized by the 
notion of an isolated horizon El [H |3 QUI]. Us ~ 
ing isolated horizons, it has been possible to derive the 
laws of black hole mechanics, use it as a basis for the 
quantum black hole entropy calculations and find unex- 
pected properties of hairy black holes in Einstein- Yang- 
Mills theory; see 1 19] and references therein. Most im- 
portantly for our purposes, isolated horizons have also 
proved to be useful in numerical relativity. For exam- 
ple, isolated horizons provide a coordinate invariant 
method of calculating the angular momentum and mass 
of a black hole |35]. They can be used to obtain bound- 
ary conditions for constructing quasi-equilibrium initial 
data sets 136LI371 13^,391. They might have a role in wave- 
form extraction [15]. A pedagogical review of isolated 
horizons from the numerical relativity perspective can 
be found in EDl . 

In this paper, we are more interested in the dynamical 
regime when the MTT is not null. A spacelike MTT con- 
sisting of future-marginally trapped surfaces is called 
a Dynamical Horizon (DH). Thus, a dynamical horizon 
is a spacelike 3-surface equipped with a given foliation 
by FMOTSs. The properties of a dynamical horizon are 
studied in detail in |9, 10, 40]. The case when the horizon 
is very close to being isolated but still evolving dynam- 
ically has been studied in |41, 42] and its Hamiltonian 
treatment is considered in |43]. Note that the local ex- 
istence of DHs follows from the local existence of MTTs 
because if ®( n ) < at any given time, it will continue 
to be strictly negative for at least a short duration. We 
elaborate on the spacelike property below. 

A timelike MTT will be called a timelike membrane 
(TLM). A TLM cannot be considered to represent the 
surface of a black hole since a time-like surface is not 
a one-way membrane, and both ingoing and outgoing 
causal curves can pass through it. In some instances, we 
shall use the term "horizon" loosely to refer to a generic 
marginal surface or a MTT without any further quali- 
fiers. The exact meaning should hopefully be clear from 
the context. 

An explicit example of a dynamical horizon is pro- 
vided by the Vaidya spacetime which describes the 
gravitational collapse of null dust |44, 45, 46]. (See 
also 1 47] for further examples in spherically symmetry). 
More generally, figure Q] depicts a dynamical horizon H 
bounded by two MOTSs Si and S2. S is a typical mem- 
ber of the foliation. The vector i a is the future directed 
unit timelike normal to H, f a is tangent to H and is 
the unit outward pointing spacelike normal to the cross- 
sections. A fiducial set of null normals is 



r = -=(t fl +f), 



n a = —(i a -r a ). (2.7) 
V2 



As before, 



10 



and ©(„) < 0. The area of a cross- 



section S will be denoted by A5 and its radius by Rg : = 
^/As/47T. A radial coordinate on H will be denoted by 
r; the cross sections of H are the constant r surfaces. The 
3-metric and extrinsic curvature of H will be denoted 




FIG. 1: A dynamical horizon H bounded by MOTSs Si and S2. 
£" is the outgoing null normal, n a is the ingoing null normal, f 
is the unit spacelike normal to the cross-sections, and i a is the 
unit timelike normal to H. E is a Cauchy surface intersecting 
H in a 2-sphere S. T" is the unit timelike normal to E and R" 
is the unit space-like outward pointing vector normal to S and 
tangent to E. 



respectively by q a i, and K fl j,, and q a {, is the 2-metric on S. 

Figure Q] shows also a Cauchy surface E intersecting 
a dynamical horizon H. This intersection S will always 
be assumed to be one of the given cross-sections of H. 
The unit timelike normal to the horizon is T" and the 
unit outward pointing spacelike normal to S within E 
is R a . The three metric and extrinsic curvature of E are 
denoted by q a ^ and K a t, respectively. The fiducial set of 
null normals to S arising naturally from E are 

£« = _L(T" + R fl ), ff = -^={T a - R a ) . (2.8) 

V2 V2 

A boost transformation of the form of equation d2.2> con- 
nects (£",n a ) and (l a ,n a ): 



l a = ft , n a = f 



-l^fl 



(2.9) 



When the horizon settles down and becomes null, an 
infinite boost (/ — ► 00) is required to go from (2 a ,n a ) to 
(£ a ,n a ). 

We conclude this sub-section with a short summary of 
some basic properties of a dynamical horizon: 

Topology: The cross-sections of a DH can be either 
spherical or toroidal 0j|l(J,|l3,|34|]. Toroidal topol- 
ogy is possible only in exceptional cases when 

V(t)ab> the scalar curvature R of S, £{&t e y R„h£ b , 
and Z, tt (defined in section HTH all vanish on S HfJl . 
We shall therefore always take the cross-sections to 
be spherical. There are no similar results for cross- 
sections of TLMs. However, we use an apparent 
horizon tracker which can only locate spherical 
AHs [48] and therefore all observed MOTSs have 
spherical topology. 

Second Law: The area of the cross-sections of a DH in- 
creases along f a i9l ll0ll . Thus, if we choose a time 
evolution vector field t" for which t ■ f > 0, then 
the area of the dynamical horizon will increase 
in time, and this result can be called the second 



law for dynamical horizons. Similarly, the area 
of a TLM decreases if ©i n \ < 0, and increases if 



®(n) > 0. 



(») 



Foliation and Uniqueness: Any given spacelike MTT 
cannot have more than one distinct dynamical 
horizon structure on it |40]. This means that a DH 
can have one, and only one foliation by FMOTSs. 
This further implies that if a Cauchy surface £ 
does not intersect a given DH in one of the pre- 
ferred cross-sections, then the intersection can- 
not be a MOTS at all. Thus, different choices of 
Cauchy surfaces in general lead to different dy- 
namical horizons. There are however some con- 
straints on the location of dynamical horizons and 
trapped surfaces as proved by Ashtekar and Gal- 
loway 1 40]. For example, they show that given a 
dynamical horizon H (along with a mild generic- 
ity assumption), there cannot be any trapped sur- 
faces (and therefore no DHs) contained entirely in 
the past domain of dependence of H. See also 
146, 491 for further discussion. 



2. The signature of a MTT 

As discussed above, MTTs have been shown to exist 
for a large and physically interesting class of MOTSs, 
and this is borne out in a large number of numeri- 
cal simulations where MOTSs are located and evolved 
smoothly. How many of these MTTs are actually dy- 
namical horizons? In other words, when is a MTT space- 
like? The first result in this direction was obtained by 
Hayward ] 16] (see also |35]). Using the Raychaudhuri 
equation for £ a , it can be shown that an MTT is space- 
like if a < 0, null if ex. = and timelike if a > 0, where 



ab 



R„d a l b 



CT {l)ab cr (i) + K ab 



(2.10) 



In writing this expression, it is assumed that t a and 
n a are extended off H geodetically, so that £„6(« is 
meaningful. The term in the numerator is strictly pos- 
itive in the case of dynamical horizons if the matter 
fields satisfy, say, the null energy condition. It vanishes 
for isolated horizons. The denominator is negative for 
the Vaidya spacetime and also for the stationary Kerr- 
Newman family. This captures the notion that as we 
go inside the black hole, the outgoing null rays become 
more and more converging. Assuming that the numer- 
ator of Eq. 12.101 is nowhere vanishing on H, the hy- 
pothesis that H is spacelike is equivalent to £„@/n < 0. 
As shown by Ben-Dov |50], this last condition is not 
satisfied for all MTTs; in Oppenheimer-Snyder collapse 
15111 . there exists a timelike world tube of FMOTSs with 
£ n & {e) > 0. 

The issue of the signature has been considered in 1 3 
There it is shown that if a MOTS S is strictly stably out 



ermost, and if the quantity cr^ ah a?t + R ah l a l b is non- 
zero somewhere on S (and assuming the null energy con- 
dition), then the MTT containing S is spacelike in a 
neighborhood of S. This result is stronger than Hay- 
ward's result (Eq. J2.10I ') and it shows clearly that the 
spacelike case is physically the most interesting because 
cr (t)ab cr 1(,\ + Rab£ a P wn l n °t vanish in a non-stationary 
situation. It also shows, somewhat surprisingly, that 
even if matter or radiation is falling into a black hole 
only in the form of say, a single narrow beam from a par- 
ticular direction, the entire MTT is spacelike. One might 
naively have thought that the MTT would be spacelike 
only on portions where the energy flux is non-zero, and 
null otherwise. This is not the case because of the elliptic 
nature of the equations governing the deformations of a 
MOTS. 3 

In all the examples we present later, it turns out that 
MOTSs form in pairs, i.e., just after a MOTS So appears 
initially, it bifurcates into "outer" and "inner" MTTs, 
H out and H,„ respectively. The initial MOTS So is the 
common cross section of H ou t and H,„, and the union 
Htot = H ou t 1J Hj n forms a single smooth manifold, as far 
as we can tell numerically (though a more detailed anal- 
ysis of the differentiability of Htot is required). In par- 
ticular, the area of the cross-sections is a differentiable 
and monotonic function on this manifold. Furthermore, 
H ou t is spacelike, even on the initial cross-section So- 
This means that the inner MTT H,„ is, by continuity, ini- 
tially spacelike in an open neighborhood of So- How- 
ever, in some cases H, n soon acquires a mixed signature 
and becomes more and more timelike, and ends up as 
a TLM. We strongly suspect that such a bifurcation is a 
general phenomenon whenever a new MOTS is formed. 
The MOTSs on the inner MTT are not strictly-stably- 
outermost and thus H,„ is not required to be spacelike 
according to the results of 1341 . 

There is one configuration where the existence of an 
inner MTT is plausible. Figure |2 shows two MOTSs 
S(i),(2) surrounded by a common MOTS S ou t', ®(£) van- 
ishes on all these surfaces. Let us assume that Sm, Sm, 
and S ou t are all strictly-stably-outermost and that de- 
forming S(i) and S(2) outward yields strictly untrapped 
surfaces S',^ and S', 2 \- Similarly, suppose that deforming 

Sout inwards gives a strictly trapped surface S' out . Then, 
since 0(f) must change sign somewhere between S' out 
and S' (11 or S', 2 \, it is plausible that there is a MOTS S,„ in 
the intermediate region inside S ou t and outside Si and 
S2. This argument is supported by a recent result by 
Schoen [52] which shows the existence of a MOTS be- 



3 The Oppenheimer-Snyder case studied in 1 50] does not satisfy the 
hypotheses of these theorems because for this case, the matter fields 
have a discontinuity at the surface of the star. Further examples in 
spherical symmetry are studied in 1 47] where the matter fields are 
smooth. 




FIG. 2: Two MOTSs Sm and S/ 2 ) surrounded by a common 
MOTS Sont- Spheres lying just inside these FMOTSs must have 
negative outgoing expansion. Thus, there must be a inner 
trapped horizon S,„ inside S ou t which encloses Sh\ and S/ 2 )- 



tween a trapped (in our case S' t ) and an untrapped sur- 
face (in our case SL, [J S| 2 A It might be possible to ex- 
tend this proof to rigorously prove the existence of Sj n 
in our case, and to check whether it is topologically a 
sphere. S(i ), S(2), and S ou t are cross sections of a dynam- 
ical horizon while S,„ is a cross-section of an MTT, not 
necessarily a dynamical horizon. 



III. APPLICATIONS 



the time coordinate t as the radial coordinate r on H by 
considering H to be embedded into the spacetime M. by 
means of the map 



F(r,6,<p) = (t = r,x i = F'{r,8,4>)). 



(3.2) 



The maps F are known as soon as the MOTSs are found 
by the AH tracker. As a frame on H we choose 



6(1) = dg, e (2 ) 



sin 9 



9<*/ 



-(3) 



= 3r. 



(3.3) 



Hence, e^) connects a point on a MOTS at a certain in- 
stant of time with a corresponding point on the MOTS 
at the next instant of time. Note that this choice of frame 
breaks down at the poles of the sphere. To apply for- 
mula <I3.H . the frame J3.3J on H must be pushed forward 
to M. by means of the embedding F in the standard way: 



e (3) = (i,a,.T 1 ,a r F 2 ,a r F 3 ). 



(3.4) 



This enables us to calculate *l(f)(/) using the 3-metric on 
the Cauchy surface, and the lapse and shift. 

Having calculated the matrix q (;)(;) anc ^ assuming 
its determinant to be positive, we can easily calculate 
the unit vector f". It is simply the outward pointing 
unit spacelike vector which is a linear combination of 
( e (i)/ e (2)/ e (3))/ and is orthogonal to em and e^)- This 
construction of f a will also work in the timelike case, but 
not in the null case where q (,)(,) becomes degenerate. 



This section discusses some possible applications of 
dynamical horizons. These ideas are illustrated using 
concrete numerical examples later in SectionHVl 



A. Computing the signature of a MTT 

From a numerical standpoint, it is more convenient 
to deduce the signature of H by directly calculating the 
induced metric q a \,, rather than from Eq. <2.10t by cal- 
culating £ n &t£\ which requires extensions of £ a and n a 
away from the horizon. The signature of H is then de- 
termined by the sign of the determinant of q a \, which 
is gauge independent; note that the determinant is it- 
self gauge dependent. To calculate q a \, we find a frame 
e?.% (i = 1,2,3) on H, i.e., three smooth vector fields on 

H which are pointwise linearly independent. We then 
simply need to compute the determinant of the matrix 



%■)(/) 



*h*h- 



(3.1) 



We construct a frame on H as follows. Let (t, x l ) 
(i = 1,2,3) be the spacetime coordinates on M. used in 
the numerical simulation. The MTT H is topologically 
I x S 2 (I some interval in R) so that we can assume co- 
ordinates (r, 6, <p) on it. Here (6, (p) are standard coor- 
dinates on S 2 and r is a radial coordinate. We can use 



B. Angular momentum and mass 

Let q> a be a rotational vector field on H tangent to 
each cross-section. 4 The angular momentum of a cross- 
section S associated with q>" is given by 

Js' ) = -^f s K a b <p a ? b d 2 V. (3.5) 

We refer to ] 10] for a justification for this formula. The 

interpretation of } s as angular momentum is most 
clear cut when q> a is a rotational symmetry on H, i.e., 
when £(pK a b = and £(pq a b = 0- See 1 35] for a method 
of finding Killing vectors suitable for numerical imple- 
mentation. Booth and Fairhurst have shown that this 
formula also arises from a Hamiltonian calculation 1 42 



As we shall see below, /g is also gauge invariant when 
q> a is only divergence free, and not necessarily a sym- 



metry vector. However, } s 
more general cp a . 



(V) 



may not be meaningful for 



4 This means that cp" is tangent to S, has closed integral curves, and 
is normalized so that its integral curves have an affine length of 2n, 
and it vanishes at exactly two points on S. 



If a cross-section S has radius R$ and angular momen- 
tum jl , we can meaningfully talk about the mass: 



M 



(V) _ 



2Rc 



n 



-4(/W)2. 



(3.6) 



This mass has the same dependence on the area and 
angular momentum as in the Kerr solution. There is a 
meaningful balance law for the mass and furthermore, it 
satisfies a physical process version of the first law I9L I10H . 
Equation 13.51 uses the metric q ab and K a f, and extrin- 
sic curvature of the dynamical horizon. It is more con- 
venient to recast this in terms of the metric q ab and ex- 
trinsic curvature K a t, of the partial Cauchy surface £ (see 
figure^. It is also more convenient to work with the null 
normals (l a , it") defined in equation J2.8I . It is clear that 
(l a ,n a ) must be related to the old null normals (£",n a ) 
by a boost transformation, i.e., there must exist a posi- 
tive function f on S such that 



oa c na 



lafl 



ft and n a =f~ 1 n 



(3.7) 



After some simple algebra, equation 13.51 can be written 
as: 



(13.91 may not look like a rotational vector field; in par- 
ticular it may vanish at more than just two points on the 
sphere even when S is close to axisymmetry 5 This is 
work in progress. The results presented below all use 
the method described in ] 35] of finding Killing vectors 
based on the Killing transport equations. This reduces 
the problem of finding Killing vectors on a sphere to 
the diagonalization of a 3 x 3 matrix, and integrating a 
1 -dimensional ordinary differential equation. We have 
found this method to be quite reliable for the cases when 
the horizon is sufficiently close to axisymmetry, even 
in cases when the coordinate system is not adapted to 
the axial symmetry. Thus, it works well for the head- 
on collision and axisymmetric neutron star collapse, but 
only at very early and late times for a non-axisymmetric 
black hole collision. This caveat only affects the exam- 
ple of section IIV Bl It is important to keep in mind that 
this Killing transport method is not reliable for check- 
ing whether the horizon is close to axisymmetry; this re- 
quires an independent calculation of £q>q a b to verify that 
it is sufficiently small. Finally, we emphasize that this 
method is also not guaranteed to produce a divergence 
free rotational vector field; this must also be checked in- 
dependently. 



r(«0 



.-^K ab R^d 2 V- 



£<p]nfd 2 V 



C. Multipole moments 



The second integral vanishes precisely when q> a is di- 
vergence free, i.e., when cp" is a symmetry of the area 
element on S. In this case: 



's 



l-^K ab R a cp b d 2 V. (3.8) 



In particular, this will be true when q> a is a symmetry 
of the metric q a },, but the divergence free condition is 
much weaker than this. For example, following 1191 . 
we can always construct a divergence free vector field 
on a 2-sphere even in the absence of axisymmetry as 
follows. Let h be any smooth function on S, and g an- 
other smooth function satisfying e djid^g = 0, where 
£ a l, is the volume form on S. It is easy to check explic- 
itly that the following vector field is automatically di- 
vergence free: 



= ge«%h. 



(3.9) 



The integral curves of q> a are the level curves of h. In par- 
ticular, if h is chosen to be a geometric quantity such as, 
say, the curvature R, and g chosen such that <p a has affine 
length In, then q> a will coincide with an axial Killing 
vector, if it exists. Therefore, q> a can be viewed as an er- 
satz axial symmetry vector field even in the absence of 
axisymmetry 

However, we haven't as yet satisfactorily imple- 
mented the above construction due to numerical diffi- 
culties arising from errors in taking derivatives of the 
scalar curvature. Furthermore, the <p a coming from eq. 



The notion of multipole moments play a very impor- 
tant role in Newtonian gravity and classical electrody- 
namics. Let us focus on classical electrodynamics in 
Minkowski space with axisymmetric charge and current 
distributions p and j a respectively, given on a sphere S of 
radius R$- Let (6, <p) be coordinates on S; p and j a , being 
axisymmetric, are functions only of 8. The electric mul- 
tipoles E n and magnetic multipoles B n are respectively 
defined as 



F - R" 



I pP„{cos6)d 2 V, 



(3.10) 



B„ 



-R 



n+l 



JxdP n {cos6)) -nd 2 V, (3.11) 



where P„ is the n Legendre polynomial, d denotes the 
standard derivative operator on a sphere, and ft is the 
unit outward normal to the sphere. For black holes, the 
analogs of the electric and magnetic multipole moments 
are respectively the mass and angular momentum mul- 
tipole moments. Motivated by this analogy, there ex- 
ist meaningful definitions of the source multipole mo- 
ments for an isolated horizon 1 53]. Roughly speaking, 
these definitions correspond to taking the moments of 
the free data on an axisymmetric isolated horizon, and 
knowledge of these moments is sufficient to construct 
the entire horizon geometry. 



5 We thank Ivan Booth for this comment. 



For dynamical horizons, we can generalize the con- 
struction of 1 53] to construct a set of multipole moments 
which capture the geometry of a dynamical horizon at 
any instant of time, and which are furthermore equal to 
the isolated horizon multipole moments when the black 
hole is isolated. The analog of charge density is (propor- 
tional to) the scalar curvature on S: 



and the angular momentum current is 



(3.12) 



(3.13) 



The moments of these quantities will give the desired 
multipole moments. We could also use q c a K cb r b instead 
of q c a K cb R above; the two expressions are related by a 
boost transformation. Just as for angular momentum, 
the final expressions for the multipole moments given 
below will be boost invariant if the <p a used in their def- 
inition is divergence free. To define the moments, we 
need a preferred coordinate system on S so that we can 
define the preferred spherical harmonics. 

The construction of the preferred coordinate system 
(6, <p) on S is the same as given in [53]: cp e [0, 2n) is 
the affine parameter along q>" and t, := cos 6 E [— 1, 1] is 
defined by the condition 



D fl £ = -j^Sbaf" ■ 
K S 



(3.14) 



The freedom to add a constant to £ is removed by requir- 
ing its integral over S to vanish: § s t,d 2 V — 0. When ap- 
plied to a Kerr black hole, these invariant coordinates 
turn out to be the same as the usual Boyer-Lindquist 
(9, <p) coordinates. 

The mass and angular multipole moments are then re- 
spectively: 



M„ 

In 



R"M S 



— I {£?„(£)} d 2 V , (3.15) 

TL JS 

^- I {e a \d b P n {Q)K ac R c )d 2 V 



^-^P' n iC)K nl ^R"d 2 V 



(3.16) 



where V' n (£) = dP„ (£) lit,. We have used equation J3~14l 
to obtain the final expression for /„ above. This form 
clarifies the relation of /„ to the angular momentum and 
also demonstrates the gauge invariance of J n when q> a is 
divergence free. Using the Gauss-Bonnet theorem, it is 
trivial to check that Mo = Mg and ]\ = Js- Jo vanishes 
because we do not consider any topological defects. Fur- 
thermore, these expressions are well suited for numeri- 
cal computation because they involve only quantities on 
the Cauchy surface and an integral over the MOTS. 



D. The energy and angular momentum fluxes 

Hawking's area theorem shows that if matter satis- 
fies the null energy condition, then the area of the event 
horizon can never decrease. This is one of the central 
results of black hole physics, and it leads to the classical 
picture of the black hole growing inexorably as it swal- 
lows matter and radiation. Therefore, one might expect 
there to be a balance law relating the increase in area to 
fluxes of matter and radiation crossing the event hori- 
zon. However, the teleological nature of event horizons 
is again a problem; there cannot exist any such local bal- 
ance law for the area of the event horizon. A clear ex- 
ample is seen in the Vaidya spacetime where the event 
horizon is formed in flat space and its area increases in 
anticipation of matter falling into the black hole at a later 
time; see 1 19] for a discussion. 

For DHs, it is possible to obtain an exact balance law 
for the area increase |9, 10]; i.e., given two cross-sections 
Si and S2 with radii R\ and R2 respectively, and with 
S2 lying to the outside of Si, the increase in the radius 
is given by the sum of the energy flux due to matter 
(F^) and gravitational radiation (J 7 ^'), both of which 
are manifestly positive.: 



R2-R1 



jr('») + jte) 



(3.17) 



where 



F {m) = I V2T ab t a £ b dRd 2 V, (3.18) 

Jh 

F {s) = ^J H {W (l) \ 2 + \tf}dRd 2 V. (3.19) 

Here \a (() \ 2 := cr {e)ab cr^ y |£| 2 := t, a t, a where £ fl is a vector 
on S defined as 



C := V2f b ? c V c £ 



b> 



(3.20) 



and d 2 V is the natural geometric volume element on H. 
The extra factors of 2 and V2 in the above equations as 
compared to the corresponding equations in 1 10], arise 
because of our normalization convention i ■ n = — 1; llOll 
uses I ■ n — —2. 

See 1 10] for additional reasons why J 7 ^ has the right 
properties to be viewed as the flux of gravitational radi- 
ation. Equation <3.17t is an exact statement about black 
holes in full non-linear general relativity, and it is the 
analog of the Bondi mass balance law at null infinity. 

From a numerical point of view, ?*•> is inconve- 
nient to calculate, especially when the horizon is settling 
down and is close to being null. First of all, we have di- 
rect access only to the fiducial null normals (2 a , n a ) de- 
fined in eq. J2.8t and not to (£", n a ) themselves. The two 
sets of null normals are related to each other by a boost 
transformation £ a = fl a , n — f~ 1 n". Under this trans- 
formation, <Tf = fcj. Similarly, it is easy to show that 



where 



r = /v 



CO 



where 



k" = fH c Vch and 



to 



= q ah n c V r l 



c*fe. 



(3.21) 



(3.22) 



Here k" and co a are tangent to the cross-sections of the 
DH. When the DH approaches equilibrium, f — » oo. 
However, the value of J-(%> itself remains finite. All 
fields with a bar remain finite even when the horizon 
becomes null even though / diverges While this is not 
a problem analytically, this does cause numerical errors 
in the transition to equilibrium when we multiply a very 
small quantity on the horizon with a very large one. This 
is consistent with the results of |41] where it is found 
that |C(i)| 2 is the most important when the horizon is 
close to equilibrium. 

Let t be the time coordinate used to label the Cauchy 
surfaces. Using this coordinate, we can identify the di- 
vergence of various terms appearing in !F^'. We start 
by rewriting !F^> as: 



fit) 



SF/JW + ltf/f*"*. (323) 



The integrand on the right hand side can be expanded 
as 



(i 



R/ 4 |k| 



Rf(\a {l) \ 2 -co-K) 



R\io\ 



(3.24) 



Let us look at the various terms in this expression. First, 
Co" can be shown to be equal to the angular momentum 
current; for an axial symmetry vector cp a , the angular 
momentum is simply the integral of (p a co a over the cross 
section of the MTT. Thus, to a need not vanish even when 
the MTT becomes an isolated horizon. The |a)| 2 term in 
the flux can, in some sense, be viewed as the flux of ro- 
tational energy entering the horizon. Now consider k" . 
For an isolated horizon, £ b V b l a cx£ a because in this case 
l a is guaranteed to be geodetic. This implies R a = 0. 
On the dynamical horizon side, we can choose suitable 
extensions of £" (and n a ) away from the MTT so that 
K a = 0. The shear Pip. on the other hand contains most 
of the non-trivial information about the radiation falling 
into the black hole. It vanishes on an isolated horizon 
as it should, and it is independent of any extensions of 
£ a , n" a way from the MTT. Therefore, in the examples of 
section HVl we shall usually plot <X/j\ to show the energy 
flux falling into the horizon. 

The angular momentum also obeys a balance law sim- 
ilar to equation i3.17\ : 



J 



(m) 



J< 



(g) 



f T ab i a tp b d 3 V, (3.26) 

Jah 

~ I P ab C f q ab d 3 V (3.27) 

167T J AH 



where P ab := K ab - Kq ab . Unlike the energy flux T^\ 
the angular momentum flux J^> is not positive definite. 
Also, j{%> vanishes when cp" is an axial Killing vector 
on H. Thus, angular momentum is conserved in the ax- 
isymmetric vacuum case, as it should be. 



IV. EXAMPLE NUMERICAL SIMULATIONS 

In this section, we apply the ideas discussed in the 
previous sections to three concrete numerical simula- 
tions: i) A head-on collision of two black holes starting 
with Brill-Lindquist initial data; ii) A non-axisymmetric 
black hole collision using puncture initial data with non- 
vanishing linear momentum and iii) Axisymmetric col- 
lapse of a neutron star. Each of these three cases is quite 
well known in the numerical relativity literature, and all 
have been well studied. This section aims to further ex- 
plore these examples using the tools described in Sec- 
tionHIH 



A. Head-on collision with Brill-Lindquist data 



The Brill-Lindquist initial data |54] for binary black 
holes represent initial data for two non-spinning black 
holes without any orbital angular momentum. The 
reader can consult a review on initial data, such as 15511 
for details. Here we simply note that these initial data 
are conformally flat and time-symmetric: 



lab 



V>%b , 



K ab = . 



(4.1) 



The manifold E is 1R 3 with two points removed (the 
punctures). The only equation to be solved is the flat 
space Laplace equation for the conformal factor: 



Aip = 0. 



(4.2) 



Let A denote the shortest distance between the two 
punctures as measured with respect to the fictitious flat 
background metric & ab ; the physical proper distance be- 
tween the punctures is actually infinite. It was shown in 
1 54] that each of the punctures is actually an asymptot- 
ically flat region. The total ADM mass of the common 
asymptotic region is 



h-h = Jo 



(m) 



■J< 



(,?) 



(3.25) 



m. 



= 2a 



(i) 



2a 



(2)' 



(4.3) 
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and the ADM masses of the two punctures are 



m 



in 



,ADM 

(i) = 


-v.. , 2a (D a (2) 
" 2a (D + d 


(4.4) 


,ADM 

(2) - 


_ o n , , 2a (D«(2) 

- 2cc {2) + d . 


(4.5) 
(4.6) 



These are exact results, irrespective of the distance d be- 
tween the punctures. In the next two sub-sections, we 
look at two different regimes (i) the far limit when d is 
large and (ii) the merger of the two holes starting from 
relatively small values of d. 



1. The far limit 

Before presenting the results from the numerical evo- 
lution of this data, it is instructive to look at a special 
case which is amenable to analytic treatment, namely, 
in the far limit where the separation between the holes 
is very large: d 3> am, Kn)- In this case, there are 
two MOTSs surrounding each of the punctures with- 
out any common MOTS surrounding them. The angu- 
lar momenta of the two black holes are trivially zero 
because the extrinsic curvature vanishes. What about 
the mass? Should m^ M and mfiy? be identified with the 
masses of the black holes? There are three difficulties 
with this. First, these ADM masses also include contri- 
butions from radiation present in the respective asymp- 
totic regions. Secondly, if this identification is correct, 
m ?iT (* = ■"•' 2 ) * s su PP ose d t° be the mass of the black 
hole for all values of d, even when the two black holes 
are very close to each other. Shouldn't the mass of the 
black holes in this regime also include, say, contribu- 
tions from the tidal distortions produced by the other 
hole? Finally, the strategy of using the asymptotic re- 
gions to define black hole masses is not applicable gen- 
erally, say in the case when there are matter fields and 
the topology of L is just R 3 , or in Misner data 1 56] where 
the two black holes do not have their own individual 
asymptotic regions. 

From the isolated /dynamical horizon perspective, 
since the black holes have zero angular momentum, 
from equation J3.8J . the irreducible mass is the correct 

measure of mass in this case: mm = \/dm /167T where 

am is the area of the MOTS around each of the punc- 
tures. Let us then calculate the mass of the black holes as 
a power series in 1/d. To simplify calculations, put the 
origin of coordinates at the location of the first puncture 
and the other puncture on the z-axis at (0,0, d). Intro- 
duce the usual spherical coordinates (r,G,(p) so that the 
conformal factor becomes explicitly 



c K r,e ) = i + a -^ + a -^(i- 2 -^ 



r2 



We see that due to axisymmetry, there is no dependence 
on <p. Let the surface of the FMOTS around the origin be 
given by the equation r = h(6). In the limit when d — > 
oo, the initial data reduces to Schwarzschild in isotropic 
coordinates so that the horizon is located at r = am. 

Higher order effects can also be explicitly calculated. 
It turns out |57] that up to 0(d~ 3 ), the location of the 
MOTS is given by 



r = 



"(1)*(2) . a (l) a (2), 
« (i) - -±L±± + -L^li(« (2) - « (1) COSf 

« (1) « (2) 



3 

5 
7 fc (D 



( a.%) ~ 3 a (i) a (2) cos 

-o(d~ 4 ) 



« m^2 (cos 



(4.8) 



where P 2 is the second Legendre polynomial. Using this 



result, the horizon mass mm = J am/16n can be calcu- 
lated and, somewhat surprisingly, the mass is the same 
as the ADM mass even up to third order: 



m (l) = 2a (l) 



?M +0 (0. (4.9) 



(4.7) 



This relation was verified numerically for a sequence of 
BL data with different values of d. However, we did not 
have sufficient resolution to estimate the leading order 
deviation between mm and m^f. Similarly, the shear 
of the horizon vanishes up to third order indicating that 
the individual horizons are isolated to an excellent ap- 
proximation. As we shall see below, the individual hori- 
zons are isolated even for relatively small values of d 
once the common MOTS has formed. 



2. Numerical results for the merger phase 

We performed a numerical evolution starting with 
Brill-Lindquist initial data. Working in units where 
the total ADM mass is unity, the punctures were lo- 
cated at z = ±0.5, and the individual black holes had 
equal masses. Thus 2am = 2oci2) — 0-5- The domain 
had an explicit octant symmetry and extended up to 
x,y,z = 96. Near the outer boundary the spatial res- 
olution was h — 1.6, and near the punctures we used 
mesh refinement to increase the resolution successively 
up to h = 0.0125, so that the individual horizon diam- 
eters contained initially 32 grid points. We used fourth 
order accurate spatial differencing operators, and a third 
order Runge-Kutta time integrator. 

We excised 1581 coordinate spheres with a radius of 
r e = 0.0625 about the punctures from the domain, cor- 
responding to a diameter of 10 grid points. We used 
the AEI BSSN formulation |58, 59] for time evolution, 
using the boundary conditions also described in 1 58]. 
These boundary conditions are known to be incompati- 
ble with the Einstein equations. We used a 1 + log slic- 
ing condition 1 60] starting from a = 1, and a zero shift. 
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FIG. 3: Coordinate shapes of the horizons at t = 1 in the xz 
plane. A common horizon has formed, and the inner and outer 
common horizons have already separated. Compare figure|5] 
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Horizon metric determinant at t=1 



This makes both the individual and the outer common 
horizon grow in coordinate space. We used the Cactus 
framework |61, 62], the Carpet mesh refinement driver 
I63ll64j l, and the CactusEinste in infrastructure. We lo- 
cated the apparent horizon surfaces with J. Thornburg's 
AHFinderDirect lisllrj^l . 

In this setup, the apparent horizon has two discon- 
nected components in the initial data, and a common 
MOTS forms shortly after t = 0.5. The individual hori- 
zons are null up to numerical errors (consistent with 
the result on the smallness of crm in the far limit), and 
their masses are essentially constant up to numerical er- 
ror. As discussed in section Hi B 2l and figure |21 the com- 
mon MOTS forms initially as a single surface but then 
bifurcates: an outer horizon which is strictly-stably- 
outermost, and an inner one which becomes strictly un- 
tapped on being deformed inwards. Figure |3] shows 
the shapes of the individual and the inner and outer 
common MOTSs at time t = 1, where the inner and 
outer common MTTs have already noticeably separated. 
As expected, the outer MTT is purely spacelike while 
the inner MTT, being spacelike initially, becomes partly 
timelike quickly. Figure||]shows the horizon world tube 
metric signature at t = 0.6 and t = 1. At later times, 
the outer MTT tends to become null (as expected), while 
the inner MTT becomes completely timelike, and then 
becomes so distorted at about t = 1.2 that it cannot be 
reliably tracked any more. This coordinate distortion is 
already evident in figure |3| and the horizon discretisa- 
tion used in the apparent horizon finder is inaccurate 
near the neck of the inner horizon 14811 . Figure |5] shows 
the time evolution of the masses M = \J As/16ti of the 
individual and the common horizons (in this case, the 
angular momentum vanishes identically). If M 00 is the 
asymptotic value of the mass of the outermost horizon 
at late times, then M ADM — M°° is, in principle, a reliable 
way of estimating the amount of energy radiated away 
to infinity in the form of gravitational waves. This dif- 
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FIG. 4: Determinant of the horizon world tube's three-metric 
vs. latitude 9 at t = 0.6 and t = 1. The individual MTTs are 
null, i.e., det q = (up to numerical errors). The common 
outer MTT is spacelike (i.e., det q > 0) and it tends to null 
at late times. The inner common MTT is partially timelike at 
t = 0.6; later it becomes completely timelike. 



ference could be used as a consistency check on other es- 
timates using the extracted waveforms at large distances 
from the black holes. However, our emphasis in this pa- 
per is on the dynamics of the merger and not on long 
duration stable evolutions. Our simulations do not last 
long enough to estimate M°° reliably. 

Another feature of the horizons, shown in figure |5j is 
that while the common outer MTT increases in area as 
expected, the area of the common inner MTT decreases 
monotonically. This is explained as follows. Initially, 
when the common MOTS is just formed, by continu- 
ity with the outer MTT, the inner MTT is spacelike for a 
very short duration (much before t = 0.6) and it is thus 
a DH for this duration. However, this DH is being tra- 
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FIG. 5: Irreducible mass vs. time for the individual and the 
common MTTs. The outer common MTT grows and accretes 
mass, while the inner MTT shrinks and loses mass. 



versed in the inwards direction (i.e., along — f a ) so that its 
area appears to decrease. Shortly after its formation, the 
inner MTT becomes partly timelike and later fully time- 
like. Recall that for a TLM, the area decreases if ©(„) < 0. 
Thus, both the spacelike and timelike portions of the in- 
ner MTT contribute to its monotonic area decrease. This 
behavior of the outer MTT is roughly similar to what 
was found in |47] for spherically symmetric horizons; 
however due to spherical symmetry, the horizons in |47] 
did not have any cross sections of mixed signature. 

Figure |6] demonstrates how the common outer ap- 
parent horizon grows. The energy flux vanishes at the 
poles, and the shear (but not the total flux) is maximum 
at the equator. The horizon is spacelike all the time, but 
it becomes increasingly isolated at late times as it ap- 
proaches equilibrium. Thus the rate of area increase be- 
comes smaller and the fluxes also becomes correspond- 
ingly smaller. 

Let us now consider the higher mass multipoles M n 
(all the ] n s vanish identically). Here, since all quanti- 
ties are symmetric with respect to a reflection about the 
equatorial plane, M n = for odd n. Figure plots the 
mass quadrupole moment M2 and also M4 and M& of 
the outer and inner common MTTs as a function of time. 
We expect that the black hole should eventually settle 
down to a Schwarzschild solution by radiating away all 
of its higher multipole moments. Clearly, for the outer 
MTT, M2, M4 and Mg all become smaller with time, ap- 
proaching zero. However, the run did not last long 
enough for us to obtain the asymptotic fall-off rate. It 
is interesting to note that, as far as we can tell, the mul- 
tipole moments for the inner MTT do not vanish asymp- 
totically This tells us that the spacetime near the inner 
MTT is not close to Schwarzschild even at late times. 
At even later times, all the inner horizons presumably 
cease to exist (see next paragraph) and the spacetime ap- 
proaches Schwarzschild everywhere. 

We conclude this section with some remarks on the 
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FIG. 6: Energy flux and shear |c7;)| 2 through the outer com- 
mon horizon vs. latitude 8 at t = 0.6, and the total energy flux 
vs. time. The shear vanishes at the poles and the black hole 
settles down exponentially. 



eventual fate of the inner MTT. First of all, as expected, 
the outer MTT eventually settles down and approaches 
future timelike infinity. The inner MTT shrinks and ap- 
proaches the two individual horizons which are essen- 
tially stationary. It is interesting to speculate on how, if 
at all, the inner MTT will merge with the two individual 
MTTs. Does the inner MTT "pinch off" into two indi- 
vidual horizons? If the inner MTT is indeed the one pre- 
dicted by [52], then it has a priori curvature bounds. If 
these curvature bounds are maintained in the limit, then 
the inner horizon cannot pinch off. It is more likely that 
the two individual MTTs merge first with each other and 
then later, perhaps also with the inner MTT. It would 
be interesting to investigate this question further. If 
the inner MTT does indeed merge smoothly with the 
two individual MTTs, then the set of all MTTs in this 
case would form one single smooth 3-manifold. Fur- 
thermore, the area of the cross-section of this manifold 
would be monotonic in the outward direction - travers- 
ing this manifold in the outward direction means going 
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FIG. 7: Some mass multipole moments vs. time for the inner 
and outer MTTs for the head-on collision. The multipole mo- 
ments for the outer horizon all approach their Schwarzschild 
values (i.e., 0) but the inner horizon does not seem to do so. 



forward in time on the individual and outer MTTs, and 
backward in time on the inner MTT. 

We are not able to settle these issues numerically in a 
conclusive manner because the inner MTT becomes so 
distorted at late times that the AH tracker is no longer 
able to track it. This is because the AH tracker can only 
locate star-shaped surfaces and, as is clear from figure 
[3J the inner MTT will not necessarily be star-shaped at 
later times. Furthermore, our gauge choice in which we 
allow the outer MTT to grow in coordinate space, makes 
the inner MTT shrink and therefore harder to resolve at 
later times. 



B. Non-axisymmetric black hole collision 

The head-on collision described above does not incor- 
porate any effects of angular momentum. In this section, 
we remove the restriction of axisymmetry by taking ini- 
tial configurations in which the black holes are orbiting 



around each other. We use the so called "puncture" data 
introduced by Brandt and Brugmann 1 66], which is a 
generalization of the Brill-Lindquist construction. The 
data is still taken to be conformally flat, but now no 
longer assumed to be time symmetric. 

We performed a numerical evolution of puncture ini- 
tial data corresponding to the innermost stable circular 
orbit as predicted in 1 67]. This model was also studied 
as "QC-0" with the Lazarus perturbative matching tech- 
nique miH and later in 0. 0. B IE Iz3l In our setup, 
the punctures were located at x = ±1.168642873, and 
their mass parameters were m = 0.453, and their mo- 
menta were py = ±0.3331917498. The domain had an 
explicit rotating quadrant symmetry and extended up 
to x,y, z = 10. Near the outer boundary the spatial res- 
olution was h = 0.4, and near the punctures we used 
mesh refinement to increase the resolution successively 
up to h — 0.025, so that the individual horizon diam- 
eters contained initially 16 grid points. We used fourth 
order accurate spatial differencing operators, and a third 
order Runge-Kutta time integrator. 

We excised |58] coordinate spheres with a radius of 
r e = 0.075 about the punctures from the domain, corre- 
sponding to a diameter of 6 grid points. We used again 
the AEI BSSN formulation |58, 59] for time evolution, a 
1 + log slicing condition 1 60] starting from a lapse that is 
one at infinity and zero at the punctures, and a T driver 
shift condition, starting from a rigid co-rotation with an 
angular velocity of to = 0.06. We also used a drift cor- 
recting shift term similar to 1 71, 72] to keep the individ- 
ual horizons centered about their initial locations. 

As previously, we used the Cactus framework 1 61, 62], 
the Carpet mesh refinement driver 1 63, 64], and the 
CactusEinstein infrastructure. We solved the initial 
data equation with M. Ansorg's TwoPuncture solver 
17311 . and we located the apparent horizon surfaces with 
J. Thornburg's AHFinderDirect |48]. 

This setup contains two initially separated horizons 
that rotate around each other for a fraction of an orbit 
before a common horizon forms |68, 7Q]. Its ADM mass 
is M ADM = 1.00788, the initial proper horizon separation 
is L ps 4.99 M ADM , and the horizons have initially the an- 
gular momentum / rs 0.78 M ADM and angular velocity 
Q f=a 0.17/M ADM . The common apparent horizon forms 
at about t = 17.5, which we verified through pretracking 
& 

Figure |H] shows the shape of the various MOTSs at a 
time f = 18, a short while after the common horizon has 
formed. The qualitative behavior of the various MTTs is 
exactly the same as in section IIV A 21 Figure El shows 
the irreducible mass of the outer and inner MTTs as a 
function of time. Again, the behavior is qualitatively the 
same as we saw in the head-on collision. 

Figure^2 snows the flux of gravitational wave energy 
falling into the outer horizon at t = 18.4 and also the 
shear \ca\ | 2 at the same time, for the outer and individ- 
ual horizons. The 2-d contour plots of the shear |<77^| 2 
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FIG. 8: Coordinate shapes of the MOTSs at t = 18 for the non- 
axisymmetric black hole collision. Note that the individual 
horizons are locked in place through the co-rotating coordi- 
nate system and through an adaptive shift condition. 



FIG. 10: A plot of the irreducible mass M, n . = V ' A/\6n as 
a function of time for the outer an inner MTTs in the non- 
axisymmetric black hole collision. As expected, the outer MTT 
has increasing area while the inner MTT shrinks. 



Horizon metric determinant at t=18 
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FIG. 9: Determinant of the MTT three-metric at t = 18. As 
in the head-on case, the outer MTT is purely spacelike while 
the inner MTT is partly spacelike and partly timelike. At later 
times, it becomes purely timelike. The individual MTTs are 
null at this time. 



and the total flux on the horizon shows in detail how 
gravitational radiation is falling into the horizon. Un- 
like in the head-on case (fig. [6j, the shear and the flux 
are now no longer axisymmetric. Therefore, the flux is 
no longer constant along the <p direction but its maxima 
still lie on the equator. The shear on the other hand, now 
has its maximum on the poles and its minima lie on the 
equator. It would be interesting to further investigate 
the behavior of \ca\ | 2 and the energy flux as a function 
of time and for different physical situations to gain a bet- 
ter understanding of how a black hole grows. 

Let us now turn to the rotational vector cp" on the 
outer horizon and the quantities such as angular mo- 
mentum, mass, and multipole moments associated with 
it. The simulation presented here was run only up to 
t fa 19 A, and the final black hole has not settled down 
sufficiently, and has not attained axisymmetry at this 
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Lie derivative of the two-metric at t=1 8 




FIG. 12: Lie derivative of the two-metric Hq,t\ a b at f = 18 on 
the <p — line. The two-metric q a f, is also shown for compar- 
ison. The quantity shown in the plots are actually the norms 

yLabiAp^ab^ and ^Zab(qab) 2 m th e coordinate system (8, <p) 
on the horizon. The vector field cp is Killing on the equator (see 
main text), but not everywhere. This shows that the horizon is 
not (yet) axisymmetric. We expect it to become axisymmetric 
at later times. Note that we have only shown the plots along 
the <p = curve and we do not have axisymmetry here. 



point. Figurell2lshows the Lie derivative of the 2-metric 
Ccp^ab on the horizon at t = 18, where cp" is the Killing 
vector candidate found by the algorithm presented in 
[35]. It is clear that Cq,q a }, is very far from at this time. 
This means that the angular momentum, mass, and mul- 
tipole moments associated with this cp" are not mean- 
ingful at this point. This is to be expected, since the fi- 
nal black hole should attain axisymmetry only on a time 
scale set by the quasi-normal mode ringdown, which 
has a period of 15.9M ADM in this case. It is interesting 
to see that our Killing vector field candidate is indeed 
Killing on the equator. This is by construction, since we 
choose the Killing vector field candidate by an integral 
along the equator; see ] 35]. However, the vector field cp" 
is far from Killing away from the equator. 

A word of caution is due here regarding the Killing 
vector finding algorithm of 1 35]. First of all, the algo- 
rithm only produces a candidate for a Killing vector, and 
an independent check is required to see whether £mf a 2, 
is sufficiently small or not. Furthermore, as mentioned 
previously, this method reduces the problem of finding a 
Killing vector on a sphere to diagonalizing a 3 x 3 matrix 
followed by integrating a 1-dimensional ODE. In partic- 
ular, the method requires that one of the eigenvalues of 
this matrix is sufficiently close to unity. While this is fine 
when the horizon is exactly axisymmetric, the subtlety 
arises when the horizon is only approximately axisym- 
metric. It is not clear how close the eigenvalue must 
be to unity for the horizon to be regarded as approxi- 
mately axisymmetric. Work is in progress to understand 
this better and to also investigate an alternate method of 
finding an appropriate cp" as discussed in section UlI Bl 



which is guaranteed to produce a divergence free vec- 
tor. 



C. Axisymmetric gravitational collapse 

1. The initial configuration 

Up to now, all of our examples have involved only 
vacuum spacetimes. In this section, we present an exam- 
ple of the gravitational collapse of a neutron star to form 
a black hole in an axisymmetric spacetime. These sim- 
ulations were performed using the Whisky code which 
deals with the matter terms of the Einstein equations in 
the framework of the Cactus toolkit. Thus, the Whisky 
code solves the conservation equations for the stress en- 
ergy tensor T ab and for the matter current density J a : 



V fl T flb = 0, V a / fl =0. 



(4.10) 



For details about the Whisky code and the implementa- 
tion of the above equations, we refer the reader to 1 75] 
and references therein. Here we shall restrict ourselves 
to describing the initial stellar configuration which is 
one of the configurations studied in 1 75]. 

The neutron star is modeled as a uniformly rotating 
ball of perfect fluid. The equation of state is taken to be a 
K = 100, r = 2 poly trope so that the pressure p and rest- 
mass density p are related according to p = Kp T . The 
equilibrium configuration is determined by the mass 
M NS , central density p c , and the angular momentum 
/ NS ; when necessary, the subscript NS is used in order to 
avoid any confusion with previously defined symbols. 
The model we take is the one denoted as "D4" in 17511 
which has M NS = 1.86M©, p c = 1.934 x 10 15 g cm" 3 , and 
/ NS = 0.543M^ S . This leads to a ratio of polar to equa- 
torial coordinate radii of 0.65, a circumferential equa- 
torial radius of 14.22 km, and a rotational frequency 
of 1295.34Hz. This equilibrium configuration turns out 
to be dynamically unstable. In practice, the instability 
is induced by uniformly reducing the pressure slightly 
throughout the star. 



2. Numerical results 

We simulated the above system on a grid with an ex- 
plicit rotating octant symmetry. The outer boundary 
was at X,y,Z — 150, and the grid spacing near the outer 
boundary was h = 3. We used mesh refinement to in- 
crease our spatial resolution in the center of the domain 
to h = 0.375 at the initial time, and progressively intro- 
duced more mesh refinement levels to increase the cen- 
tral resolution up to h = 0.046875 as the neutron star 
collapsed, based on the maximum density in the star 
1 76, 77]. We also apply third order Kreiss-Oliger dissipa- 
tion 1 78] to the spacetime (but not the hydrodynamics) 
variables. 
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FIG. 13: The average coordinate radius and the area radius as 
a function of time for the outer and inner MTTs for the neutron 
star collapse. The inner horizon is not to be trusted after t w 
140 due to lack of resolution, since its coordinate radius has 
become very small by that time. 



We find an apparent horizon starting at about t = 130; 
this time is mainly dependent on the details of how 
the collapse is induced and has no intrinsic meaning. 
The horizon is born with an irreducible mass of about 
Mj rr = 1.51 and an angular momentum of / = 0.89 
(a = 0.38), giving it a total mass of Mjf = 1.54. Some 
time after t = 185, a singularity forms in the spacetime, 
and the simulation aborts because we do not use exci- 
sion inside the apparent horizon. As before, a pair of 
MOTSs is formed, an outer and an inner one. The outer 
MTT is spacelike, has increasing area, and tends to null 
at late times. In this case, the inner MTT remains space- 
like. However, its area decreases because we are travers- 
ing it in the inward direction; in other words, the time 
evolution vector t " is such that at the inner MTT, t • f < 
so that that the area decreases along t". Our gauge con- 
ditions are such that the outer horizon grows in coordi- 
nate space while the inner horizon shrinks. After about 
t = 140, the inner horizon is so small that we do not 
have enough resolution to track it beyond that time. See 
figure ^] The areal radius of the outer MTT increases 
but not as rapidly as the coordinate radius; it levels off 
at later times. The area radius of the inner horizon de- 
creases initially and shows an increase at later times, but 
this is probably just a numerical artefact due to poor res- 
olution at later times. 

Figure ^] shows the determinant of the metric on the 
MTTs. The outer MTT is initially spacelike, which is con- 
sistent with its growing, and exponentially approaches 
null at late times. After about t = 160, the simula- 
tion cannot distinguish the horizon world tube signa- 
ture from null any more. As an example we also show 
the determinant as a function of the latitude 6 at t sa 138, 
and the average value of the determinant over the hori- 
zon as a function of time. The inner MTT is also space- 
like and becomes more and more null at least as long as 
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FIG. 14: Average of the determinant of the horizon world 
rube's three metric vs. time, and vs. latitude 9 at t = 138.24 
for the inner and outer horizons for the neutron star collapse. 



we are able to track it reliably. 

Figure ^] shows the outer horizon has grown at t = 
155 to an irreducible mass of Mj rr = 1.80 and an angu- 
lar momentum of / = 1.93 (a = 0.55), giving it a total 
mass of Mh = 1.87. For comparison, the correspond- 
ing ADM quantities are Madm = 1-86 and / ADM = 1.88 
(a = 0.54). Because the spacetime is axially symmet- 
ric, gravitational waves cannot carry away angular mo- 
mentum. That means that the spin a = ] / M 2 is ap- 
proximately correct at late times. Unlike in the non- 
axisymmetric black hole collision discussed earlier, the 
present case is explicitly axisymmetric and there are no 
problems with locating the rotational symmetry vector. 

Figure ^] shows the mass quadrupole moment Mi 
and the angular momentum octopole moment of the 
outer and inner MTTs as a function of time. Given that 
we know the asymptotic values of the area and angu- 
lar momentum of these MTTs (the ADM values), we can 
also calculate the expected values of Mi and J$ at late 
times. The plots clearly show that the values of Mi and 
Js approach the Kerr values at later times (though this 
matching is not exact, presumably due to numerical er- 
rors). Also note that Mi is noisy. We have observed such 
noise only in simulations that include matter, and we 
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FIG. 15: The total mass Mn„ and angular momentum / as a 
function of time for the outer and inner MTTs for the neutron 
star collapse. 



find that this noise is much improved by using artificial 
dissipation on the spacetime variables (which we do). 
The angular momentum multipoles seem unaffected. 



V. DISCUSSION 

In this article, we have applied the dynamical hori- 
zon formalism to numerical simulations of black hole 
spacetimes. The main theme in this formalism is to 
take trapped surfaces seriously as a way of describing 
black hole physics. Marginally trapped surfaces behave 
more regularly that one might have expected previously 
and they are useful for extracting interesting physical 
information about the horizon. We have shown how 
the mass, angular momentum, multipole moments, and 
the flux of energy due to in-falling gravitational radi- 
ation and matter can be calculated in a coordinate in- 
dependent way (given a particular time slicing of our 
spacetime). We have implemented these ideas numeri- 
cally and shown three concrete examples. In these ex- 
amples, we see how the black hole is formed, how it 
grows, and how it settles down to an isolated Kerr black 
hole. We have also seen that the dynamical horizon for- 



FIG. 16: Horizon mass quadrupole moment M2 and angular 
momentum octopole moment J$ vs. time for the neutron star 
collapse. For comparison, the values for a Kerr black hole with 
the same ADM mass and angular momentum as the initial 
data are also shown. 



malism is valuable for exploring the geometry of the 
trapped region. It allows us to classify various types of 
trapped surfaces which might appear during the course 
of a gravitational collapse or a black hole coalescence. 
Finally, these ideas can also be viewed as a set of di- 
agnostic tools which allow us to keep track of what is 
going on during the course of a numerical simulation, 
and whether numerical results make sense and satisfy 
some basic, but non-trivial properties in the strong field 
region. 

Some suggestions for future work: 

i. As mentioned in the text, the calculation of the ax- 
ial vector cp" for non-axisymmetric cases is not yet 
satisfactory. We have used the method suggested 
in 13511 which works well enough at early and late 
times, when the horizon is approximately axisym- 
metric. However, in general, the result is not guar- 
anteed to be divergence free and thus the angular 
momentum not guaranteed to be gauge invariant. 
We have not yet implemented the generalization 
described in section fTll Bl satisfactorily: this is work 
in progress. 
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ii. The accuracy of the numerical examples that we have 
shown decreases with time, and this is a common 
feature of most present day black hole numerical 
simulations. Thus, we have not been conclusively 
able to prove that the black hole settles down to 
Kerr (though there are strong indications that this 
does happen). We have not been able to extract the 
rate at which equilibrium is reached, thereby ex- 
tending Price's law (see |79] and e.g. 1 80]) to more 
general situations, but this is, in principle possible 
and requires more stable and accurate simulations. 
Similarly, we have not been able to accurately cal- 
culate the asymptotic value of the black hole mass 
M°°. The difference M ADM - M°° is, in principle, 
a reliable estimate of the amount of energy radi- 
ated to infinity. While the ADM mass is hard to 
calculate reliably during the simulation because of 
the finite grid and low resolution in the asymp- 
totic region, it can usually be calculated accurately 
from the initial data themself . Calculating M°° and 
understanding this estimate of the radiated energy 
requires more accurate and stable runs, applied to 
diverse and realistic initial data. The results of |41] 
could also be used to study the approach to equi- 
librium. 

iii. It would be useful numerically to have a gauge con- 
dition which ensures that the horizon stays at the 
same coordinate location at all times. While such 
conditions are not difficult to find in the isolated 
case, dynamical situations are harder. Given the 
location of an outer MOTS at a particular instant of 
time, the results and methods of 1 34] can be used 
to predict the location of the MOTS at the next in- 
stant by solving an elliptic equation on the MOTS. 
This could be used to construct appropriate gauge 
conditions and evolution schemes which take the 
horizon geometry into account 1 49, 81]. 

iv. What happens to the inner horizon of figures |3J 
andlSf As described in section ll V A 21 the eventual 
fate of these inner MTTs and the two individual 
horizons is not yet known, and would be inter- 
esting to investigate further. This requires simu- 
lations with higher resolution near the inner hori- 
zons, different gauge conditions, and perhaps also 



AH trackers capable of handling non-star-shaped 
surfaces, and perhaps also higher genus surfaces. 

v. Can the methods of |34] be extended for MOTSs 
which are not strictly-stably-outermost? In this re- 
gard, it would be interesting to study the stabil- 
ity operator L^ introduced in [34]. For a strictly- 
stably-outermost MOTS, the principle eigenvalue 
of Lj] turns out to be strictly positive and this is 
an important ingredient in the existence results. 
A numerical computation of the eigenvalues of 
this operator, especially during the transition be- 
tween inner and outer MTTs and for the inner non- 
spacelike MTTs might lead to further insights. 
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